Accessibility

This section will introduce calculating accessibility metrics in R. Specifically, we’ll work to understand questions such as:

• What is the travel time between specific origin and destination points?
• Which locations are within 30 minutes of a health facility?
• How can we measure changes in accessibility?

Data sources

• OpenStreetMap: Public available information on road networks, buildings, etc. In short, a Wikipedia for maps.
• Google Maps and Mapbox provide APIs for querying travel time between origin and destinations and by travel time (driving, walking, public transit); incorporates live traffic information.
• Other private sector companies like TomTom and INRIX provide more granular information (eg, distribution of speeds across road segments, both in real time and histroic); however, data sources can be more expensive to access.

Today, we’ll focus approaches using OpenStreetMaps. We’ll primarily rely on Open Source Routing Machine (OSRM) for determining accessibility, which is a routing engine that uses OpenStreetMap data.

Setup

Copy/paste the following code into a new RStudio script, replacing “YOURFOLDERPATHHERE” with the folder within which you’ll place this R project:

library(usethis)

use_course(
  url = "https://github.com/dime-worldbank/gis-transport-training/archive/main.zip",
  destdir = "YOURFOLDERPATHHERE"
)

Setup

Install new packages

install.packages(c("sf",
                   "leaflet",
                   "geosphere"),
                 dependencies = TRUE)

And load them

library(here)
library(tidyverse)
library(units)     # Define units
library(sf)        # Simple features
library(leaflet)   # Interactive map
library(osmdata)   # Query OSM data
library(osrm)      # Routing using OpenStreetMaps

Example travel time between two schools

Load school data and convert to spatial object

schools_df <- read.csv(here("DataWork", 
                            "DataSets", 
                            "Final",
                            "schools.csv"))

schools_sf <- st_as_sf(schools_df,
                       coords = c("longitude", "latitude"),
                       crs = 4326)

Example travel time between two schools

Map two schools

leaflet() %>%
  addTiles() %>%
  addCircles(data = schools_sf[1:2,], opacity = 1, weight = 10)

Example travel time between two schools

Determine driving time between schools (measured in minutes)

tt_drive <- osrmTable(src = schools_sf[1,],
                      dst = schools_sf[2,])

tt_walk <- osrmTable(src = schools_sf[1,],
                     dst = schools_sf[2,],
                     osrm.profile = "foot")

print(tt_drive$durations)

    2
1 2.5

print(tt_walk$durations)

   2
1 20

Example travel time between two schools

Determine optimal route

drive_sf <- osrmRoute(src = schools_sf[1,],
                      dst = schools_sf[2,])

walk_sf <- osrmRoute(src = schools_sf[1,],
                     dst = schools_sf[2,],
                     osrm.profile = "foot")

leaflet() %>%
  addTiles() %>%
  addPolylines(data = drive_sf, opacity = 1, weight = 10) %>%
  addPolylines(data = walk_sf, color = "orange", opacity = 1)

Example travel time between two schools

Exercise:

1. Determine the driving and walking travel time between two different schools
2. Determine geometry of the shortest driving and walking routes
3. Add the travel time as a variable on the geometry objects
4. Map the shortest path geometries

Example travel time between two schools

Solution:

1. Determine geometry of the shortest driving and walking routes
2. Determine the driving and walking travel time between two different schools and add the travel time as a variable on the geometry objects

## School ids
id_1 <- 10
id_2 <- 20

## Route geometries
drive_sf <- osrmRoute(src = schools_sf[id_1,],
                      dst = schools_sf[id_2,])

walk_sf <- osrmRoute(src = schools_sf[id_1,],
                     dst = schools_sf[id_2,],
                     osrm.profile = "foot")

## Add travel time as a variable
drive_sf$travel_time <- osrmTable(src = schools_sf[id_1,],
                                  dst = schools_sf[id_2,])$durations

walk_sf$travel_time <- osrmTable(src = schools_sf[id_1,],
                                 dst = schools_sf[id_2,],
                                 osrm.profile = "foot")$durations

print(c(drive_sf$travel_time, walk_sf$travel_time))

[1]  5.9 41.2

Example travel time between two schools

Solution:

3. Map the shortest path geometries

## Map
leaflet() %>%
  addTiles() %>%
  addPolylines(data = drive_sf, opacity = 1, weight = 10,
               popup = ~paste0("Travel time: ", travel_time, " minutes")) %>%
  addPolylines(data = walk_sf, color = "orange", opacity = 1,
               popup = ~paste0("Travel time: ", travel_time, " minutes"))

Travel time between multiple origins and destinations

Compute travel time from the first school to the next 99 schools

tt_drive <- osrmTable(src = schools_sf[1,],
                      dst = schools_sf[2:99,])

hist(tt_drive$durations)

Travel time between multiple origins and destinations

Compute optimal routes:

1. osrmRoute only allows 1 O-D pair, so we apply over destinations
2. osrmRoute is slow, so we just look at 10 O-D pairs

drive_sf <- map_df(2:11, function(i){
  osrmRoute(src = schools_sf[1,],
            dst = schools_sf[i,])
})

Travel time between multiple origins and destinations

leaflet() %>%
  addTiles() %>%
  addPolylines(data = drive_sf) %>%
  addCircles(data = schools_sf[2:11,], color = "red", opacity = 1) %>%
  addCircles(data = schools_sf[1,], color = "green", weight = 10, opacity = 1)

Travel time between multiple origins and destinations

Exercise: A student lives near school 1, but has options to go to all other schools. The student will have to walk to school and doesn’t want to cross a major road (trunk or motorway) to get to the school. How many schools can the student consider?

(To minimize computational time, only consider the first 20 schools).

Hint: Load the road files and subset to major roads. Which routes don’t intersect with a major road?

roads_sf <- read_sf(here("DataWork", 
                        "DataSets", 
                        "Final",
                        "roads.geojson"))

roads_sf <- roads_sf %>%
  filter(highway %in% c("trunk", "motorway"))

# Make one row
roads_sf <- roads_sf %>% st_union()

Travel time between multiple origins and destinations

Solution:

walk_sf <- map_df(2:20, function(i){
  osrmRoute(src = schools_sf[1,],
            dst = schools_sf[i,],
            osrm.profile = "foot")
})

walk_sf$inter_mjr_road <- st_intersects(walk_sf, roads_sf, sparse = F)

table(walk_sf$inter_mjr_road)

TRUE 
  19 

Isocrhones

iso_sf <- osrmIsochrone(
  loc = schools_sf[1,],
  breaks = seq(from = 0, to = 20, length.out = 5)
)

head(iso_sf)

Simple feature collection with 4 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 36.62947 ymin: -1.364402 xmax: 36.93593 ymax: -1.13673
Geodetic CRS:  WGS 84
  id isomin isomax                       geometry
1  1      0      5 MULTIPOLYGON (((36.784 -1.2...
2  2      5     10 MULTIPOLYGON (((36.84123 -1...
3  3     10     15 MULTIPOLYGON (((36.84123 -1...
4  4     15     20 MULTIPOLYGON (((36.81645 -1...

Isocrhones

leaflet() %>%
  addTiles() %>%
  addCircles(data = schools_sf[1,], color = "red") %>%
  addPolygons(data = iso_sf)

Isocrhones

Exercise:

1. Make an Isochrone from the centroid of the city. Maximum distance should be 40 minutes
2. Map the Isochrone and schools
3. Challenge What proportion of schools are within 10 minutes of the center of the city? (Hint: Use st_join to perform a spatial join).

city_sf <- read_sf(here("DataWork", 
                        "DataSets", 
                        "Final",
                        "city.geojson"))
head(city_sf)

Simple feature collection with 6 features and 13 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 36.67803 ymin: -1.370704 xmax: 36.99025 ymax: -1.234921
Geodetic CRS:  WGS 84
# A tibble: 6 × 14
  GID_2   GID_0 COUNTRY GID_1 NAME_1 NL_NAME_1 NAME_2 VARNAME_2 NL_NAME_2 TYPE_2
  <chr>   <chr> <chr>   <chr> <chr>  <chr>     <chr>  <chr>     <chr>     <chr> 
1 KEN.30… KEN   Kenya   KEN.… Nairo… <NA>      Dagor… <NA>      <NA>      Const…
2 KEN.30… KEN   Kenya   KEN.… Nairo… <NA>      Dagor… <NA>      <NA>      Const…
3 KEN.30… KEN   Kenya   KEN.… Nairo… <NA>      Embak… <NA>      <NA>      Const…
4 KEN.30… KEN   Kenya   KEN.… Nairo… <NA>      Embak… <NA>      <NA>      Const…
5 KEN.30… KEN   Kenya   KEN.… Nairo… <NA>      Embak… <NA>      <NA>      Const…
6 KEN.30… KEN   Kenya   KEN.… Nairo… <NA>      Embak… <NA>      <NA>      Const…
# ℹ 4 more variables: ENGTYPE_2 <chr>, CC_2 <chr>, HASC_2 <chr>,
#   geometry <POLYGON [°]>

center_sf <- city_sf %>% st_union() %>% st_centroid()

Isocrhones

Solution:

## Make IsoChrone
iso_sf <- osrmIsochrone(
  loc = center_sf,
  breaks = c(0, 10, 20, 30, 40)
)

## Map
leaflet() %>%
  addTiles() %>%
  addPolygons(data = iso_sf) %>%
  addCircles(data = schools_sf, color = "red")

Isocrhones

Solution:

## Spatial join
schools_iso_sf <- schools_sf %>%
  st_join(iso_sf)

## Proportion within 10 minutes
mean(schools_iso_sf$isomax <= 10, na.rm = T)

[1] 0.05945337

Isocrhones

Instead of using Isochrones, we could also also determine the percent of schools within 10 minutes using osrmTable. One disadvatange of this approach is that the code can take longer to run; a separate query to OSRM is made for each school.

osrmTable(src = center_sf,
          dst = schools_sf)

Isodistance

Isochrones show polygons that are reachable within a certain time. We can also compute an Isodistance, which shows polygons that are reachable within a certain driving distance.

iso_sf <- osrmIsodistance(
  loc = schools_sf[1,],
  breaks = c(2000, 4000, 6000, 8000, 10000) # meters
)

palette_1 <- colorNumeric(
  palette = "viridis", domain = iso_sf$isomin
)

leaflet() %>%
  addTiles() %>%
  addCircles(data = schools_sf[1,], color = "red", opacity = 1) %>%
  addPolygons(data = iso_sf,
              fillColor = ~palette_1(isomin),
              fillOpacity = 0.7,
              stroke = F)

Isocrhones

The osrmIsochrone function only allows determining accessibility from one distance. In many cases were are interested in understanding accessibility to the closest location among many locations; for example, travel time to the nearest school.

Strategy:

1. Make grid across city
2. Compute fastest travel time for each grid to each school
3. Make map of travel times across grids

Isocrhones

1. Make grid across city

hex_sf <- st_make_grid(
  x = city_sf,
  cellsize = 0.03, # decimal degrees, approximately 5km
  square = F # when F, makes hexagonal grid
) %>%
  st_as_sf()

hex_sf <- st_intersection(hex_sf, 
                          city_sf %>% st_union())

leaflet() %>%
  addTiles() %>%
  addPolygons(data = hex_sf)

Isocrhones

2. Compute fastest travel time for each grid to each school. To minimize computation time, we pick 5 random schools.

schools_lim_sf <- schools_sf %>%
  arrange(runif(n())) %>%
  head(5)

hex_sf$tt_nearest_school <- lapply(1:nrow(hex_sf), function(i){

  hex_sf_i <- hex_sf[i,] %>% st_centroid()
  
  min_tt <- osrmTable(src = hex_sf_i, 
                      dst = schools_lim_sf)$durations %>%
    min()
  
  return(min_tt)
}) %>%
  unlist()

Isocrhones

3. Make map of travel times across grids

palette_1 <- colorNumeric(
  palette = "viridis", domain = hex_sf$tt_nearest_school
)

leaflet() %>%
  addTiles() %>%
  addPolygons(data = hex_sf,
              fillColor = ~palette_1(tt_nearest_school),
              fillOpacity = 0.7,
              stroke = F) %>%
    addCircles(data = schools_lim_sf, color = "red", opacity = 1) 

Routing considering live traffic

Google:

• The Distance Matrix API allows querying travel times considering: (1) live traffic conditions, (2) typical traffic conditions for a date/time (e.g., typical travel time on a Monday at 3pm), and (3) different travel modes (driving, walking, biking, transit).
• The gmapsdistance package facilitates using the API. The gmapsdistance function works similarly to the osrmTable function:

library(gmapsdistance)
results <- gmapsdistance(origin = c(lat, lon), 
                        destination = c(lat, lon), 
                        mode = "driving",
                        key = "GOOGLE-API-KEY") 

Routing considering live traffic

Mapbox:

• The Mapbox Directions API works similarly to Google’s Distance Matrix API. It allows querying travel times considering (1) live traffic conditions, (2) typical traffic conditions for a date/time, and (3) different travel modes (driving, walking, biking).
• The mapboxapi package facilitates using the API. The mb_matrix function works similarly to osrmTable function:

library(mapboxapi)
results <- mb_matrix(origins = point1_sf, 
                        destinations = point2_sf, 
                        profile = "driving-traffic",
                        access_token = "MAPBOX-API-KEY") 

Routing considering alternate scenerios

Routing using a graph

1. Obtaining a road network file (e.g., shapefile of roads)
2. Manipulating the network (eg, changing speeds, removing a road)
3. Converting the road network to a graph (eg, with nodes and weighted edges) and routing along the graph.

Resources for this approach:

• tidygraph: R package
• GOSTnets: Python package (developed by World Bank)

Routing considering alternate scenerios

Routing using a gridded approach

One disadvantage with the above approach is that converting the road network spatial file into a graph can be cumbersome. Another approach that can be more straightforward but assumes some simplifcation of the road network:

1. Obtaining a road network file (e.g., shapefile of roads)
2. Manipulating the network (eg, changing speeds, removing a road)
3. Create a grid across the network (e.g., 100 meter grids).
4. For each grid cell, determine the time it takes to cross the grid using the speed limit of the fastest road that intersects with the grid.
5. For grid cells that don’t interesect with a road, assuming some walking speed
6. Use the grid for routing between grid cells

Resources for this approach:

• gdistance: R package

Thank you